library(tidyverse)
library(ggplot2)
library(grid)

plot_path        <- "output/plots"

main <- function() {
  
  df <- expand.grid(X = seq(0, 10, .1), Y = seq(0, 10, .1))
  df$Z <- df$X * df$Y / 100

  ggplot(df, aes(x = X, y = Y, z = Z)) + 
    geom_contour_filled(show.legend = FALSE) +
    annotate(geom = "rect", xmin = 0, xmax = 0.5, ymin = 0, ymax = 9.975, color = "black", fill = "black") +
    annotate(geom = "rect", xmin = 0, xmax = 9.975, ymin = 0, ymax = 0.5, color = "black", fill = "black") +
    labs(x = "Distance from causally correct\nspecification at the estimand", 
         y = "Distance between\ngiven causal summary\nand closest estimable\ncausal summary") +
    coord_cartesian(expand = FALSE, ylim = c(0, 11.5), xlim = c(0, 12)) +
    scale_x_continuous(breaks = c(0)) +
    theme_bw() +
    theme( 
      legend.direction      = "vertical",
      legend.box.spacing    = unit(0.1, "pt"),
      legend.title.align    = 0.5,
      legend.box.background = element_rect(colour = "black", size = 0.8),
      legend.title          = element_text(size = 14, face = "bold"),
      legend.text           = element_text(size = 12),
      axis.title            = element_text(size = 14),
      axis.title.y = element_text(angle = 0, vjust = 0.5),
      axis.text.y = element_blank(),
      axis.text.x = element_text(size = 14, hjust = 1.5),
      axis.text             = element_blank(),
      axis.ticks.x          = element_blank(),
      axis.ticks.y          = element_blank(),
      panel.grid.major      = element_blank(),
      panel.grid.minor      = element_blank(),
      panel.background      = element_blank(),
      axis.line             = element_line(colour = "black"),
      panel.border          = element_blank(),
      axis.line.y = element_line(arrow = grid::arrow(length = unit(0.3, "cm"), 
                                                      ends = "last")),
      axis.line.x = element_line(arrow = grid::arrow(length = unit(0.3, "cm"), 
                                                     ends = "last"))) +
    scale_fill_viridis_d(option = "magma", begin = 0.2, end = 1.0) +  
    geom_point(aes(color = Z), alpha = 0) +
    scale_color_viridis_c(option = "magma", name= "Bound on\nthe error", 
                          begin = 0.0, end = 1.0, breaks = c(0.01, 0.45, 0.85),
                          labels = c("Zero", "Small", "Large"))
  
  ggsave(sprintf("%s/isobias.pdf", plot_path), width = 8, height = 6, units = "in")

}

main()
